%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% File:               log_likelihood3.m
%
% Authors:            Sergio Ascencio and Miguel Rueda
%
% Description:        Log likelihood function for two step ML estimator (3 players).
%
% Language:           MATLAB R2013b (8.2.0.701) 64 Bit
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [LL,g,H,auxLL]=log_likelihood3(beta,Y_PRI,Y_PAN,Y_PRD,X_PRI,X_PAN,X_PRD)

betas=reshape(beta,[size(X_PAN,2),6]);

%Base outcome is L (minimum representation) order of coefficients L, M, H
betas_PRI=[zeros(size(X_PRI,2),1) betas(:,1:2)];
betas_PAN=[zeros(size(X_PAN,2),1) betas(:,3:4)];
betas_PRD=[zeros(size(X_PAN,2),1) betas(:,5:6)];

y_PRD=dummyvar(Y_PRD+1);
y_PAN=dummyvar(Y_PAN+1);
y_PRI=dummyvar(Y_PRI+1);

p_PRD=exp(X_PRD*betas_PRD)./kron(ones(1,size(y_PRD,2)),sum(exp(X_PRD*betas_PRD),2));
p_PAN=exp(X_PAN*betas_PAN)./kron(ones(1,size(y_PAN,2)),sum(exp(X_PAN*betas_PAN),2));
p_PRI=exp(X_PRI*betas_PRI)./kron(ones(1,size(y_PRI,2)),sum(exp(X_PRI*betas_PRI),2));

LL=-sum(sum(y_PAN.*log(p_PAN)+y_PRI.*log(p_PRI)+y_PRD.*log(p_PRD),2));
auxLL=-sum(sum(y_PAN.*log(p_PAN)+y_PRI.*log(p_PRI)+y_PRD.*log(p_PRD),2).^2);

%Gradient
auxPRD_M=[-p_PRD(:,2) (1-p_PRD(:,2)) -p_PRD(:,2)];
auxPAN_M=[-p_PAN(:,2) (1-p_PAN(:,2)) -p_PAN(:,2)];
auxPRI_M=[-p_PRI(:,2) (1-p_PRI(:,2)) -p_PRI(:,2)];
auxPRD_H=[-p_PRD(:,3) -p_PRD(:,3) (1-p_PRD(:,3))];
auxPAN_H=[-p_PAN(:,3) -p_PAN(:,3) (1-p_PAN(:,3))];
auxPRI_H=[-p_PRI(:,3) -p_PRI(:,3) (1-p_PRI(:,3))];

g=-sum([X_PRI.*kron(ones(1,size(X_PRI,2)),sum(y_PRI.*auxPRI_M,2)) X_PRI.*kron(ones(1,size(X_PRI,2)),sum(y_PRI.*auxPRI_H,2)) ...
    X_PAN.*kron(ones(1,size(X_PAN,2)),sum(y_PAN.*auxPAN_M,2)) X_PAN.*kron(ones(1,size(X_PAN,2)),sum(y_PAN.*auxPAN_H,2)) ...
    X_PRD.*kron(ones(1,size(X_PRD,2)),sum(y_PRD.*auxPRD_M,2)) X_PRD.*kron(ones(1,size(X_PRD,2)),sum(y_PRD.*auxPRD_H,2))],1);

%Hessian
auxPRI_M2=X_PRI.*kron(ones(1,size(X_PRI,2)),p_PRI(:,2));
auxPRI_1M2=X_PRI.*kron(ones(1,size(X_PRI,2)),(1-p_PRI(:,2)));
auxPRI_H2=X_PRI.*kron(ones(1,size(X_PRI,2)),p_PRI(:,3));
auxPRI_1H2=X_PRI.*kron(ones(1,size(X_PRI,2)),(1-p_PRI(:,3)));

auxPAN_M2=X_PAN.*kron(ones(1,size(X_PAN,2)),p_PAN(:,2));
auxPAN_1M2=X_PAN.*kron(ones(1,size(X_PAN,2)),(1-p_PAN(:,2)));
auxPAN_H2=X_PAN.*kron(ones(1,size(X_PAN,2)),p_PAN(:,3));
auxPAN_1H2=X_PAN.*kron(ones(1,size(X_PAN,2)),(1-p_PAN(:,3)));

auxPRD_M2=X_PRD.*kron(ones(1,size(X_PRD,2)),p_PRD(:,2));
auxPRD_1M2=X_PRD.*kron(ones(1,size(X_PRD,2)),(1-p_PRD(:,2)));
auxPRD_H2=X_PRD.*kron(ones(1,size(X_PRD,2)),p_PRD(:,3));
auxPRD_1H2=X_PRD.*kron(ones(1,size(X_PRD,2)),(1-p_PRD(:,3)));

H_PRI=[-auxPRI_M2'*auxPRI_1M2 auxPRI_M2'*auxPRI_H2;auxPRI_M2'*auxPRI_H2 -auxPRI_H2'*auxPRI_1H2];
H_PAN=[-auxPAN_M2'*auxPAN_1M2 auxPAN_M2'*auxPAN_H2;auxPAN_M2'*auxPAN_H2 -auxPAN_H2'*auxPAN_1H2];
H_PRD=[-auxPRD_M2'*auxPRD_1M2 auxPRD_M2'*auxPRD_H2;auxPRD_M2'*auxPRD_H2 -auxPRD_H2'*auxPRD_1H2];


H=-[H_PRI zeros(size(H_PRI,1),size([H_PAN H_PRD],2));zeros(size(H_PAN,1),size(H_PRI,2)) H_PAN zeros(size(H_PAN,1),size(H_PRD,2));zeros(size(H_PRD,1),size([H_PRI H_PAN],2)) H_PRD];
end

